#ifndef GF2M_MATRIX_H
#define GF2M_MATRIX_H

#include <cstring>
#include <exception>
#include <memory>

namespace GF2M
{
    class GF2M_MATRIX_EXCEPTION : public std::exception
    {
    private:
        std::string reason;
    public:
        GF2M_MATRIX_EXCEPTION()
        {
            reason = "GF2M Matrix Error: Reason not specified.";
        }

        explicit GF2M_MATRIX_EXCEPTION(const std::string& s)
        {
            reason = "GF2M Matrix Error: " + s;
        }

        const char * what() const noexcept override
        {
            return reason.c_str();
        }

    };

    enum class GF2M_MATRIX_TYPE{
        GF2M_MATRIX_ZERO,
        GF2M_MATRIX_UNIT,
        GF2M_MATRIX_UTRIANGLE,
        GF2M_MATRIX_LTRIANGLE
    };

    template<typename T>
    class GF2M_MATRIX {
    private:
        unsigned int Rows;
        unsigned int Cols;
        T * Matrix;
        T ** Matrix_Row_Ptr;
        static const size_t GF2M_MATRIX_BYTES{sizeof(T)};
        static const size_t GF2M_MATRIX_BITS{sizeof(T) * 8};
    public:
        GF2M_MATRIX()
        {
            Rows = 1;
            Cols = 1;

            Matrix = new T[1];
            *Matrix = 0;

            Matrix_Row_Ptr = new T*[1];
            Matrix_Row_Ptr[0] = Matrix;
        }

        GF2M_MATRIX(const GF2M_MATRIX<T> &M)
        {
            unsigned int pseudoRows, pseudoCols;

            Rows = M.Rows;
            Cols = M.Cols;
            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
            pseudoRows = (Rows%GF2M_MATRIX_BITS==0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;

            Matrix = new T[pseudoRows * pseudoCols];
            Matrix_Row_Ptr = new T*[pseudoRows];
            std::memcpy(Matrix, M.Matrix, pseudoRows * pseudoCols * GF2M_MATRIX_BYTES);
            for(unsigned int i = 0; i < pseudoRows; i++) {
                Matrix_Row_Ptr[i] = M.Matrix_Row_Ptr[i] + (Matrix - M.Matrix);
            }
        }

        GF2M_MATRIX(unsigned int rows, unsigned int cols)
        {
            unsigned int pseudoRows, pseudoCols;

            Rows = rows == 0 ? 1 : rows;
            Cols = cols == 0 ? 1 : cols;

            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
            pseudoRows = (Rows%GF2M_MATRIX_BITS==0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;

            Matrix = new T[pseudoRows * pseudoCols];
            std::memset(Matrix, 0, pseudoRows * pseudoCols * GF2M_MATRIX_BYTES);

            Matrix_Row_Ptr = new T*[pseudoRows];
            for(unsigned int i = 0; i < pseudoRows; i++) {
                Matrix_Row_Ptr[i] = Matrix + i * pseudoCols;
            }
        }

        GF2M_MATRIX(unsigned int rows, unsigned int cols, GF2M_MATRIX_TYPE type)
        {
            unsigned int pseudoRows, pseudoCols;
            unsigned int i, j, m, n;
            T tmp;

            Rows = rows == 0 ? 1 : rows;
            Cols = cols == 0 ? 1 : cols;

            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
            pseudoRows = (Rows%GF2M_MATRIX_BITS==0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;

            Matrix = new T[pseudoRows * pseudoCols];
            std::memset(Matrix, 0, pseudoRows * pseudoCols * GF2M_MATRIX_BYTES);
            Matrix_Row_Ptr = new T*[pseudoRows];
            for(i = 0; i < pseudoRows; i++) {
                Matrix_Row_Ptr[i] = Matrix + i * pseudoCols;
            }

            switch(type) {
                case GF2M_MATRIX_TYPE::GF2M_MATRIX_UNIT:
                    for(i = 0; i < Rows; i++) {
                        m = i / GF2M_MATRIX_BITS;
                        n = i - m * GF2M_MATRIX_BITS;
                        for(j = 0; j < pseudoCols; j++) {
                            if(m == j) {
                                tmp = 1;
                                ((T *)Matrix_Row_Ptr[i])[m] = tmp<<n;
                            }
                        }
                    }
                    break;
                case GF2M_MATRIX_TYPE::GF2M_MATRIX_UTRIANGLE:
                    for(i = 0; i < Rows; i++) {
                        m = i / GF2M_MATRIX_BITS;
                        n = i - m * GF2M_MATRIX_BITS;
                        for(j = 0; j <pseudoCols; j++) {
                            if(j == m) {
                                tmp = ~0;
                                ((T *)Matrix_Row_Ptr[i])[m] = tmp<<n;
                            }
                            else if(j > m) {
                                tmp = ~0;
                                ((T *)Matrix_Row_Ptr[i])[j] = tmp;
                            }
                        }
                    }
                    break;
                case GF2M_MATRIX_TYPE::GF2M_MATRIX_LTRIANGLE:
                    for(i = 0; i < Rows; i++) {
                        m = i / GF2M_MATRIX_BITS;
                        n = i - m * GF2M_MATRIX_BITS;
                        for(j = 0; j <pseudoCols; j++) {
                            if(j == m) {
                                tmp = ~0;
                                ((T *)Matrix_Row_Ptr[i])[m] = tmp>>(GF2M_MATRIX_BITS-n-1);
                            }
                            else if(j < m) {
                                tmp = ~0;
                                ((T *)Matrix_Row_Ptr[i])[j] = tmp;
                            }
                        }
                    }
                    break;
                default:
                    break;
            }
        }

        ~GF2M_MATRIX()
        {
            delete[] Matrix_Row_Ptr;
            delete[] Matrix;
        }

        unsigned int resize(unsigned int rows, unsigned int cols) {
            unsigned int pseudoRows, pseudoCols;
            if(rows == 0 || cols == 0) {
                return -1;
            }
            delete[] Matrix_Row_Ptr;
            delete[] Matrix;
            Matrix_Row_Ptr = nullptr;
            Matrix = nullptr;
            Rows = rows;
            Cols = cols;
            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
            pseudoRows = (Rows%GF2M_MATRIX_BITS==0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;

            Matrix = new T[pseudoRows * pseudoCols];
            std::memset(Matrix, 0, pseudoRows * pseudoCols * GF2M_MATRIX_BYTES);

            Matrix_Row_Ptr = new T*[pseudoRows];
            for(unsigned int i = 0; i < pseudoRows; i++) {
                Matrix_Row_Ptr[i] = Matrix + i * pseudoCols;
            }
        }

        inline unsigned int getPseudoRows() const
        {
            return (Rows%GF2M_MATRIX_BITS == 0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;
        }

        inline unsigned int getPseudoCols() const
        {
            return (Cols%GF2M_MATRIX_BITS == 0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
        }

        inline unsigned int getRows() const
        {
            return Rows;
        }

        inline unsigned int getCols() const
        {
            return Cols;
        }

        inline void set(unsigned int row, unsigned int col)
        {
            unsigned int m, n;
            if (row >= Rows || col >= Cols) {
                return;
            }
            m = col / GF2M_MATRIX_BITS;
            n = col - m * GF2M_MATRIX_BITS;
            Matrix_Row_Ptr[row][m] = Matrix_Row_Ptr[row][m] | ((T)1) << n;
        }

        inline void clear(unsigned int row, unsigned int col)
        {
            unsigned int m, n;
            if (row >= Rows || col >= Cols) {
                return;
            }
            m = col / GF2M_MATRIX_BITS;
            n = col - m * GF2M_MATRIX_BITS;
            Matrix_Row_Ptr[row][m] = Matrix_Row_Ptr[row][m] & (~(((T)1) << n));
        }

        inline int get(unsigned int row, unsigned int col) const
        {
            unsigned int m, n;
            if (row >= Rows || col >= Cols) {
                return -1;
            }
            m = col / GF2M_MATRIX_BITS;
            n = col - m * GF2M_MATRIX_BITS;
            return (Matrix_Row_Ptr[row][m] & (((T)1) << n)) == 0 ? 0 : 1;
        }

        bool isZero() const
        {
            unsigned i, j, m, n, pseudoCols;
            T tmp;
            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
            m = (Cols-1) / GF2M_MATRIX_BITS;
            n = (Cols-1) - m * GF2M_MATRIX_BITS;
            for(i = 0; i < Rows; i++) {
                for(j = 0; j < pseudoCols; j++) {
                    if(j < m) {
                        if(Matrix_Row_Ptr[i][j] != 0) {
                            return false;
                        }
                    }
                    else {
                        tmp = ~0;
                        if((Matrix_Row_Ptr[i][j] & (tmp>>(GF2M_MATRIX_BITS-n-1))) != 0) {
                            return false;
                        }
                    }
                }
            }
            return true;
        }

        bool isUnit() const
        {
            unsigned i, j, m, n, pseudoCols;
            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);

            for(i = 0; i < Rows; i++) {
                m = i / GF2M_MATRIX_BITS;
                n = i - m * GF2M_MATRIX_BITS;
                for(j = 0; j < pseudoCols; j++) {
                    if((m != j) && (Matrix_Row_Ptr[i][j] != 0)
                    || (m == j ) && (Matrix_Row_Ptr[i][m] != ((T)1<<n))) {
                        return false;
                    }
                }
            }
            return true;
        }

        GF2M_MATRIX<T> & Trans()
        {
            unsigned int rowBlocks, colBlocks;
            unsigned int pseudoRows, pseudoCols;
            unsigned int i, j, k, l;
            T BUFF;

            auto TMP = new GF2M_MATRIX<T>(Cols, Rows);

            rowBlocks = (Rows%GF2M_MATRIX_BITS==0) ? Rows/GF2M_MATRIX_BITS : ((Rows/GF2M_MATRIX_BITS)+1);
            colBlocks = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);

            for(i = 0; i < rowBlocks; i++) {
                for(j = 0; j <colBlocks; j++) {
                    for(k = 0; k < GF2M_MATRIX_BITS; k++) {
                        for (l = 0; l < GF2M_MATRIX_BITS; l++) {
                            if ((this->Matrix_Row_Ptr[i * GF2M_MATRIX_BITS + l][j] & ((T)1<<k)) != 0) {
                                TMP->Matrix_Row_Ptr[j*GF2M_MATRIX_BITS+k][i]
                                        = TMP->Matrix_Row_Ptr[j*GF2M_MATRIX_BITS+k][i] ^ ((T) 1 << l);
                            }
                        }
                    }
                }
            }

            delete[] Matrix_Row_Ptr;
            delete[] Matrix;
            Matrix_Row_Ptr = nullptr;
            Matrix = nullptr;
            Rows = TMP->Rows;
            Cols = TMP->Cols;
            pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
            pseudoRows = (Rows%GF2M_MATRIX_BITS==0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;

            Matrix = new T[pseudoRows * pseudoCols];
            Matrix_Row_Ptr = new T*[pseudoRows];
            std::memcpy(Matrix, TMP->Matrix, pseudoRows * pseudoCols * GF2M_MATRIX_BYTES);
            for(i = 0; i < pseudoRows; i++) {
                Matrix_Row_Ptr[i] = TMP->Matrix_Row_Ptr[i] + (Matrix - TMP->Matrix);
            }

            delete TMP;

            return *this;
        }

        GF2M_MATRIX<T> & operator=(const GF2M_MATRIX<T>& M)
        {
            unsigned int pseudoRows, pseudoCols;
            if(this != &M) {
                delete[] Matrix_Row_Ptr;
                delete[] Matrix;
                Matrix_Row_Ptr = nullptr;
                Matrix = nullptr;
                Rows = M.Rows;
                Cols = M.Cols;
                pseudoCols = (Cols%GF2M_MATRIX_BITS==0) ? (Cols/GF2M_MATRIX_BITS) : ((Cols/GF2M_MATRIX_BITS)+1);
                pseudoRows = (Rows%GF2M_MATRIX_BITS==0) ? Rows : ((Rows/GF2M_MATRIX_BITS)+1)*GF2M_MATRIX_BITS;

                Matrix = new T[pseudoRows * pseudoCols];
                Matrix_Row_Ptr = new T*[pseudoRows];
                std::memcpy(Matrix, M.Matrix, pseudoRows * pseudoCols * GF2M_MATRIX_BYTES);
                for(unsigned int i = 0; i < pseudoRows; i++) {
                    Matrix_Row_Ptr[i] = M.Matrix_Row_Ptr[i] + (Matrix - M.Matrix);
                }
            }
            return *this;
        }

        friend GF2M_MATRIX<T> operator+(const GF2M_MATRIX<T> &A, const GF2M_MATRIX<T> &B)
        {
            unsigned int pseudoCols;

            if(A.Rows != B.Rows || A.Cols != B.Cols) {
                throw GF2M_MATRIX_EXCEPTION("Different Matrix Size.");
            }

            GF2M_MATRIX<T> Sum(A.Rows, A.Cols);

            pseudoCols = A.getPseudoCols();

            for(unsigned int i = 0; i < Sum.Rows; i++) {
                for(unsigned int j = 0; j < pseudoCols; j++) {
                    Sum.Matrix_Row_Ptr[i][j]
                            = A.Matrix_Row_Ptr[i][j] ^ B.Matrix_Row_Ptr[i][j];
                }
            }

            return Sum;
        }

        friend GF2M_MATRIX<T> operator-(const GF2M_MATRIX<T> &A, const GF2M_MATRIX<T> &B)
        {
            return operator+(A, B);
        }

        friend GF2M_MATRIX<T> operator+(const GF2M_MATRIX<T> &A)
        {
            return A;
        }

        friend GF2M_MATRIX<T> operator-(const GF2M_MATRIX<T> &A)
        {
            return A;
        }

        friend GF2M_MATRIX<T> operator*(const GF2M_MATRIX<T> &A, const GF2M_MATRIX<T> &B)
        {
            unsigned int i, j, k, m, n;
            unsigned int pseudoCols;

            if(A.Cols != B.Rows) {
                throw GF2M_MATRIX_EXCEPTION("Matrix size not matching for multiplication operation.");
            }

            pseudoCols = B.getPseudoCols();

            GF2M_MATRIX<T> Mul(A.Rows, B.Cols);

            for(i = 0; i < A.Rows; i++) {
                for(j = 0; j < A.Cols; j++) {
                    m = j / GF2M_MATRIX_BITS;
                    n = j - m * GF2M_MATRIX_BITS;
                    if((A.Matrix_Row_Ptr[i][m] & ((T)1) << n) != 0) {
                        for(k = 0; k < pseudoCols; k++) {
                            Mul.Matrix_Row_Ptr[i][k]
                                    = Mul.Matrix_Row_Ptr[i][k] ^ B.Matrix_Row_Ptr[j][k];
                        }
                    }
                }
            }

            return Mul;
        }

        friend GF2M_MATRIX<T> operator/(const int &d, const GF2M_MATRIX<T> &M)
        {
            unsigned int pseudoCols;
            unsigned int i, j, k, m, n;
            unsigned int row, row_one;
            T * ptmp;

            if(M.Rows != M.Cols) {
                throw GF2M_MATRIX_EXCEPTION("No inverse matrix available.");
            }

            GF2M_MATRIX<T> Inv(M.Rows, M.Cols);
            if(d == 0) {
                return Inv;
            }

            GF2M_MATRIX<T> TMP = M;

            pseudoCols = TMP.getPseudoCols();

            for(i = 0; i < Inv.Rows; i++) {
                Inv.Matrix_Row_Ptr[i] = Inv.Matrix + i * pseudoCols;
                m = i / GF2M_MATRIX_BITS;
                n = i - m * GF2M_MATRIX_BITS;
                Inv.Matrix_Row_Ptr[i][m]
                        = Inv.Matrix_Row_Ptr[i][m] | ((T)1) << n;
            }

            row = 0;
            for(i = 0; i < TMP.Cols; i++) {
                m = i / GF2M_MATRIX_BITS;
                n = i - m * GF2M_MATRIX_BITS;
                row_one = TMP.Cols;
                for(j = row; j < TMP.Rows; j++) {
                    if((TMP.Matrix_Row_Ptr[j][m] & (((T)1)<<n)) != 0) {
                        row_one = j;
                        break;
                    }
                }
                if(row_one != TMP.Cols) {
                    if(row_one != row) {
                        ptmp = TMP.Matrix_Row_Ptr[row];
                        TMP.Matrix_Row_Ptr[row] = TMP.Matrix_Row_Ptr[row_one];
                        TMP.Matrix_Row_Ptr[row_one] = ptmp;
                        ptmp = Inv.Matrix_Row_Ptr[row];
                        Inv.Matrix_Row_Ptr[row] = Inv.Matrix_Row_Ptr[row_one];
                        Inv.Matrix_Row_Ptr[row_one] = ptmp;
                    }
                    for(j = row + 1; j < TMP.Rows; j++) {
                        if((TMP.Matrix_Row_Ptr[j][m] & (((T)1) << n)) != 0) {
                            for(k = m; k < pseudoCols; k++) {
                                TMP.Matrix_Row_Ptr[j][k]
                                        = TMP.Matrix_Row_Ptr[j][k] ^ TMP.Matrix_Row_Ptr[row][k];
                            }
                            for(k = 0; k < pseudoCols; k++) {
                                Inv.Matrix_Row_Ptr[j][k]
                                        = Inv.Matrix_Row_Ptr[j][k] ^ Inv.Matrix_Row_Ptr[row][k];
                            }
                        }
                    }
                    row = row + 1;
                }
                else {
                    throw GF2M_MATRIX_EXCEPTION("No inverse matrix available.");
                }
            }

            i = TMP.Rows-1;
            while(i > 0) {
                m = i / GF2M_MATRIX_BITS;
                n = i - m * GF2M_MATRIX_BITS;
                for(j = 0; j < i; j++) {
                    if ((TMP.Matrix_Row_Ptr[j][m] & (((T)1) << n)) != 0) {
                        for(k = m; k < pseudoCols; k++) {
                            TMP.Matrix_Row_Ptr[j][k]
                                    = TMP.Matrix_Row_Ptr[j][k] ^ TMP.Matrix_Row_Ptr[i][k];
                        }
                        for(k = 0; k < pseudoCols; k++) {
                            Inv.Matrix_Row_Ptr[j][k]
                                    = Inv.Matrix_Row_Ptr[j][k] ^ Inv.Matrix_Row_Ptr[i][k];
                        }
                    }
                }
                i--;
            }

            return Inv;
        }

        friend GF2M_MATRIX<T> operator/(const GF2M_MATRIX<T> &A, const GF2M_MATRIX<T> &B)
        {
            return A*(1/B);
        }

    };
}

#endif //GF2M_MATRIX_H
